1. https://github.com/dlab-geo/r-geospatial-workshop

• Click Clone or Download and download the zip file
• Upzip the zip file and make a note of the folder in which it is located

1. Open RStudio and start a new script

2. Follow along by opening r-geospatial-workshop-pt1-feb2018.html in a browser window

3. Install required libraries in RStudio

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

• Brief Intro to geospatial data

• Classes & methods for working with spatial data in R

• Mapping geospatial data

• Processing geospatial data

• Practice

Patty Frontiera

dlab.berkeley.edu

Who are you?

Why are you here?

1. Get the workshop files

https://github.com/dlab-geo/r-geospatial-workshop

• Click Clone or Download and download the zip file

• Upzip the zip file and make a note of the folder in which it is located

1. Open RStudio and start a new script

2. Follow along by opening r-geospatial-workshop-pt1-feb2018.html

1. In RStudio, install required libraries in RStudio

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

are data about locations on or near the surface of the Earth.

convey geographic information but not specific location

represent location more specifically with coordinates

46.130479, -117.134167

Coordinates only make sense when associated with a CRS!

Define:

• the shape of the Earth

• the origin

• the units

• the relationship between the system and the real world

Define: 1. the shape: the major & minor axes of an ellipsoid 2. the origin: the equator and prime meridian 3. the datum: the fit of the CRS to the Earth 4. units: eg lat/lon expressed as decimal degrees or DMS

Because of variations in 1-3, there are many geographic CRSs!

The World Geodetic System of 1984 is the default geographic CRS used today.

WGS84 is the default CRS for most GIS software

• Almost all data with lon/lat coordinates are assumed to be WGS84 unless otherwise specified

NAD83 is another common geographic CRS used by US agencies like the Census.

WGS84 and NAD83 are so similar that differences are often ignored except for applications requiring high locational accuracy.

Why do you need to know?

A Projected CRS applies a map projection to a Geographic CRS

Map Projection: mathematial transformation from curved to flat surface

All introduce distortion,

• in shape, area, distance, direction, or combo

• the larger the area the greater the distortion

No one map projection best for all purposes

Selection depends on location, extent and purpose

Spatial data is a more generic term that is not just for geographic data.

• dynamically determine spatial metrics like area and length,

• characteristics like distance and direction,

• relationships like inside and intersects from these data, and

• link data by spatial location!

Vector and Raster Data

Points, lines and Polygons

Regular grid of cells (or pixels)

• We won't be covering Raster data in this workshop*

• Why special software?

Most software can only store, edit, query and analyze numbers and text

• Databases, spreadsheets, statistical software packages

software that can import, create, store, edit, visualize and analyze geospatial data

• as geometric data types referenced to the surface of the earth via CRSs

• and methods to operate on those representations

We call software for working with geospatial data GIS

Geographic Information System

This term is commonly associated with desktop software applications.

Desktop GIS - ArcGIS, QGIS

Software extended to support geospatial data - Tableau, R

Spatial Databases - Postgresql

Web GIS - ArcGIS Online, CARTO

Custom software - leaflet web maps

You already use R

Reproducibility

Free & Open Source

Cutting edge

Thousands of Cool add-ons

• Shiny, rleaflet

Vector points, lines & polygons:

• CSV
• ESRI Shapefile
• GeoJSON
• KML

Raster grids

• TIFF, GeoTiff
• JPEG, JPEG2000
• ASCII Grid

See GIS file formats

This is one of the most, if not the most common spatial vector data file formats.

Old but everywhere!

Gotchas: 2GB limit, 10 char column names, poor unicode support

The simplest, most common format for point data

"ID","name”,”x”,”y”,”taste","price","crowded","food"
"1","babette",-122.255374,37.868428,10,4,1,1
"2","musical",-122.260698,37.868383,7,3.25,1,1
"3","starbucks",-122.266057,37.870441,6,2.95,1,0
"4","yalis",-122.266385,37.873528,7,2.95,0,0
"5","berkeleyesp",-122.268681,37.873664,3,3.25,1,0
"6","fertile",-122.268863,37.874934,5,3.25,0,1

• What columns contain the geographic data?
• What is the CRS for these data?

"ID","name”,”x”,”y”,”taste","price","crowded","food"
"1","babette",-122.255374,37.868428,10,4,1,1
"2","musical",-122.260698,37.868383,7,3.25,1,1
"3","starbucks",-122.266057,37.870441,6,2.95,1,0
"4","yalis",-122.266385,37.873528,7,2.95,0,0
"5","berkeleyesp",-122.268681,37.873664,3,3.25,1,0
"6","fertile",-122.268863,37.874934,5,3.25,0,1

• What's different?
• Does it matter?

"ID","name”,”x”,”y”,”taste","price","crowded","food"
"1","babette",122.255374,37.868428,10,4,1,1
"2","musical",122.260698,37.868383,7,3.25,1,1
"3","starbucks",122.266057,37.870441,6,2.95,1,0
"4","yalis",122.266385,37.873528,7,2.95,0,0
"5","berkeleyesp",122.268681,37.873664,3,3.25,1,0
"6","fertile",122.268863,37.874934,5,3.25,0,1

• What columns contain the geographic data?
• What is the CRS for these data?

There are many approaches to and packages for working with geospatial data in R.

One approach is to keep it simple and store geospatial data in a data frame.

This approach is most common when

• the data are point data in CSV files and

• you want to map rather than spatially manipulate or analyze the data

San Francisco Open Data Portal https://data.sfgov.org

SF Property Tax Rolls

This data set includes the Office of the Assessor-Recorder’s secured property tax roll spanning from 2007 to 2016.

We are using this as a proxy for home values.

Set your working directory first to the folder where you downloaded the workshop files!

sfhomes <- read.csv('data/sf_properties_25ksample.csv') 
head(sfhomes,6)

Set your working directory first to the folder where you downloaded the workshop files!

sfhomes <- read.csv('data/sf_properties_25ksample.csv')  
str(sfhomes)

## 'data.frame':    25000 obs. of  19 variables:
##  $ FiscalYear      : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ SalesDate       : Factor w/ 5389 levels "1920-01-30","1920-01-31",..: 4135 NA 5081 NA 3092 4626 NA 755 NA NA ...
##  $ Address         : Factor w/ 24696 levels "0000 0000                       0000",..: 3679 4927 11353 16997 7808 18856 8581 20144 11321 22546 ...
##  $ YearBuilt       : int  1923 1909 1915 1947 1981 1916 1923 2003 1925 1890 ...
##  $ NumBedrooms     : int  0 3 2 0 3 5 4 3 4 6 ...
##  $ NumBathrooms    : int  0 3 2 0 3 5 4 3 4 6 ...
##  $ NumRooms        : int  7 10 6 5 6 9 8 7 10 11 ...
##  $ NumStories      : int  1 2 1 1 2 2 2 4 2 3 ...
##  $ NumUnits        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ AreaSquareFeet  : int  2250 5830 1350 1162 1463 2304 4440 1860 3330 2246 ...
##  $ ImprovementValue: int  412556 360604 558386 344050 133307 328453 665081 314062 107655 116024 ...
##  $ LandValue       : int  618835 414696 558386 516080 189442 766392 297525 314062 41295 71386 ...
##  $ Neighborhood    : Factor w/ 40 levels "","Bayview Hunters Point",..: 39 31 11 35 26 35 39 28 11 19 ...
##  $ Location        : Factor w/ 20416 levels "(37.7081695281752, -122.449582481739)",..: 2749 17122 15070 9269 17261 7014 5152 2161 15041 11347 ...
##  $ SupeDistrict    : int  7 2 1 4 1 4 7 10 1 9 ...
##  $ totvalue        : int  1031391 775300 1116772 860130 322749 1094845 962606 628124 148950 187410 ...
##  $ SalesYear       : int  2010 NA 2014 NA 1995 2012 NA 1950 NA NA ...
##  $ lat             : num  37.7 37.8 37.8 37.8 37.8 ...
##  $ lon             : num  -122 -122 -122 -123 -122 ...

class(sfhomes)   # what is the data object type?
dim(sfhomes)     # how many rows and columns
str(sfhomes)     # display the structure of the object
head(sfhomes)    # take a look at the first 10 records
summary(sfhomes) # explore the range of values
summary(sfhomes$totvalue) # explore the range of values for one column
hist(sfhomes$totvalue)    # histogram for the totvalue column

Questions:

• What columns make these data geospatial?
• Are they vector or raster data?
• What type of geometry? Points, lines, polygons, grid cells?
• What is the CRS of these data?

Use the R base plot function to create a simple map

plot(sfhomes$lon, sfhomes$lat) # using base plot function

Use the R base plot function to create a simple map

plot(sfhomes$lon, sfhomes$lat) # using base plot function

Most widely used plotting library in R

Not specifically for geospatial data

But can be used to make fabulous maps

Great choice if you already know ggplot2

Load the library

library(ggplot2)

What's different from the plot() map?

library(ggplot2)

ggplot() + geom_point(data=sfhomes, aes(lon,lat)) + coord_map()

What's different from the plot() map?

ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1) + coord_map()

Allows you to associate a map projection with geographic coord data.

ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1) #+ coord_map()

ggplot() + geom_point(data=sfhomes, aes(lon,lat, col=totvalue)) + 
  coord_map()

What's happening here?

sfhomes_low2high <- sfhomes[order(sfhomes$totvalue, decreasing = FALSE),]

ggplot() + 
  geom_point(data=sfhomes_low2high, aes(lon,lat, col=totvalue)) + 
  coord_map()

Try it - Does the output map look different from previous one?

The order of the data in the data frame changes the map display!

Map the sfhomes data in decreasing order by totvalue.

sfhomes_high2low <- sfhomes[order(sfhomes$totvalue, decreasing = T),]
ggplot() + geom_point(data=sfhomes_high2low, aes(lon,lat, col=totvalue)) + 
  coord_map()

What does this code do?

sfhomes2010_15 <- subset(sfhomes_low2high, as.numeric(SalesYear) > 2009)

ggplot() +
  geom_point(aes(lon, lat, col=totvalue), data = sfhomes2010_15 )  +
  facet_wrap(~ SalesYear)

Visual spatial analysis!

It's a good idea to explore data values before and during mapping. Sometimes maps aren't the best way or only way to visualize data.

ggplot(sfhomes2010_15, aes(totvalue)) + 
  geom_histogram() + 
  facet_wrap(~SalesYear)


boxplot(sfhomes$totvalue~sfhomes$SalesYear)

Create basemaps on which you can display your data.

Geocode place names and addresses to get point coordinates.

and more…

Load the libary

library(ggmap)

## Google Maps API Terms of Service: http://developers.google.com/maps/terms.

## Please cite ggmap if you use it: see citation("ggmap") for details.

Some functionality may require you to register a Google API key See this StackOverflow discussion

#devtools::install_github("dkahle/ggmap")
#library(ggmap)
#register_google(key="XXXXX") # your key here

Fetch map data (default=Google) to plot using get_map

Use ?get_map for details

sf_map <- get_map("San Francisco, CA")  

This use of the command leverages the Google Geocoding API to determine the coordinates of the named location.

Fetch map data (default=Google) to plot using get_map

Use ?get_map for details

sf_map <- get_map("San Francisco, CA")  

## Source : https://maps.googleapis.com/maps/api/staticmap?center=San+Francisco,+CA&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&key=AIzaSyDHsYNZmqmTNvFn36tq470iRyXgQEU-PSE

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=San%20Francisco%2C%20CA&key=AIzaSyDHsYNZmqmTNvFn36tq470iRyXgQEU-PSE

ggmap(sf_map)

ggmap(sf_map)

Syntax similar to ggplot

# ggplot() +

ggmap(sf_map) + 
  geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))

ggmap(sf_map) +
  geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))

Create a basemap zoomed to the extent of our data

# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes$lon), lat = mean(sfhomes$lat))
sf_ctr  # take a look

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

Get a basemap zoomed to the extent of the 2015 data

# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes15$lon), lat = mean(sfhomes15$lat))
sf_ctr  # take a look

##        lon        lat 
## -122.43289   37.75964

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=37.759644,-122.432893&zoom=12&size=640x640&scale=1&maptype=terrain&language=en-EN&key=AIzaSyDHsYNZmqmTNvFn36tq470iRyXgQEU-PSE

Create a ggmap that uses customized sf_basemap

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

Finally, let's add another geospatial data layer to our ggmap.

You can use this method to add as many layers as you want to a ggmap or ggplot.

Use the read.csv function to read in a file of Bart Station locations.

bart <- read.csv("./data/bart.csv")
landmarks <- read.csv("./data/landmarks.csv")
# head (bart)
# head(landmarks)

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red") +
  geom_point(data=landmarks, aes(x=X,y=Y), col="green")

## Warning: Removed 34 rows containing missing values (geom_point).

## Warning: Removed 4 rows containing missing values (geom_point).

Questons?

There are limits to what you can do with geospatial data stored in a data frame.

The ESRI Shapefile is the most common file format for geospatial data.

Can't read in or directly plot data in a shapefile with ggplot or ggmap.

These data are not in geographic coordinates - longitude and latitude.

You can't map these with ggmap.

Spatial data objects and methods are needed to answer questions like

• What properties are in the Noe Valley neighborhood?

• What is the average property value in each SF neighborhood?

• What is the area of each SF Neighborhood and the property density?

• What properties are within walking distance (.25 miles) of the Mission neighborhood?

Classes and Methods for Spatial Data

The SP package is most commonly used to construct and manipulate spatial data objects in R.

Hundreds of other R packages that do things with spatial data typically build on SP objects.

Load the library sp

Take a look at the different types of spatial object classes supported by sp

library(sp)
getClass("Spatial") 

## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

We use the S*DF objects most frequently!

Let's transform the sfhomes15 data frame to an sp object of type SpatialPointsDataFrame

Use the sp::coordinates() method

• Sets or retrieves spatial coordinates

When transforming a DF to SPDF you are setting the coordinates.

Note - the transformation happens in place!

sfhomes15_sp <- sfhomes15  # Make a copy - why?

class(sfhomes15_sp)  # check the class of the object

coordinates(sfhomes15_sp) <- c('lon','lat')  #  Make it spatial - ORDER MATTERS!!

class(sfhomes15_sp)  # check the class of the object

sfhomes15_sp <- sfhomes15  # Make a copy - why?

class(sfhomes15_sp)  # check the class of the object

## [1] "data.frame"

coordinates(sfhomes15_sp) <- c('lon','lat')  #  Make it spatial - ORDER MATTERS!!

class(sfhomes15_sp)  # check the class of the object

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

You transformed a data frame to an SPDF using:

coordinates(sfhomes15_sp) <- c('lon','lat')

Now try this:

coordinates(sfhomes15_sp)

coordinates(sfhomes15_sp)

##             lon      lat
## 24    -122.4741 37.73601
## 35    -122.3965 37.75920
## 76    -122.3942 37.77181
## 79    -122.4754 37.71518
## 85    -122.4210 37.75876
## 90    -122.4149 37.79362
## 92    -122.4370 37.74020
## 98    -122.4437 37.78881
## 223   -122.4754 37.71492
## 229   -122.3919 37.78868
## 273   -122.3968 37.77945
## 341   -122.5040 37.75515
## 363   -122.4041 37.76422
## 366   -122.3934 37.77764
## 388   -122.4166 37.74931
## 451   -122.4112 37.73683
## 467   -122.4507 37.76448
## 473   -122.3990 37.75450
## 487   -122.4195 37.76649
## 534   -122.4225 37.78676
## 548   -122.4410 37.70854
## 549   -122.5009 37.74485
## 550   -122.4220 37.73048
## 559   -122.4947 37.75832
## 567   -122.3942 37.77181
## 595   -122.4187 37.80469
## 624   -122.4019 37.77323
## 631   -122.4229 37.80586
## 646   -122.4225 37.78676
## 658   -122.3919 37.78868
## 665   -122.4261 37.74349
## 668   -122.3904 37.78139
## 683   -122.4458 37.74398
## 774   -122.3883 37.73754
## 776   -122.4767 37.75347
## 864   -122.4326 37.74229
## 974   -122.4380 37.72624
## 1018  -122.4452 37.76906
## 1097  -122.4937 37.77999
## 1155  -122.5047 37.75390
## 1164  -122.4498 37.75451
## 1167  -122.4055 37.80424
## 1188  -122.4194 37.75596
## 1242  -122.4088 37.77953
## 1273  -122.4331 37.76019
## 1374  -122.3918 37.77899
## 1379  -122.4070 37.75674
## 1487  -122.4236 37.77800
## 1598  -122.4758 37.75497
## 1600  -122.4419 37.74582
## 1618  -122.5030 37.74587
## 1668  -122.3989 37.78990
## 1698  -122.4236 37.77800
## 1747  -122.3972 37.76314
## 1755  -122.4434 37.75159
## 1768  -122.4242 37.80551
## 1804  -122.4549 37.76402
## 1828  -122.4320 37.78195
## 1833  -122.4405 37.71825
## 1913  -122.4146 37.74714
## 1999  -122.3942 37.77181
## 2024  -122.4676 37.71048
## 2033  -122.4249 37.76548
## 2066  -122.4243 37.75206
## 2073  -122.4097 37.77613
## 2096  -122.3996 37.73859
## 2248  -122.4045 37.79929
## 2258  -122.4442 37.78322
## 2285  -122.3919 37.78868
## 2314  -122.4026 37.77879
## 2329  -122.4368 37.78228
## 2339  -122.4240 37.74366
## 2355  -122.4126 37.79035
## 2416  -122.4579 37.72632
## 2428  -122.4089 37.79286
## 2436  -122.4409 37.80451
## 2443  -122.4127 37.78995
## 2500  -122.4357 37.76468
## 2508  -122.4666 37.76347
## 2537  -122.4870 37.74605
## 2538  -122.3942 37.77181
## 2627  -122.4219 37.77389
## 2692  -122.3958 37.77705
## 2717  -122.4116 37.77806
## 2734  -122.4290 37.72977
## 2793  -122.4764 37.73390
## 2914  -122.4183 37.76712
## 2969  -122.4711 37.72241
## 2975  -122.4393 37.78018
## 3127  -122.4393 37.78951
## 3157  -122.3702 37.72838
## 3205  -122.4398 37.78746
## 3262  -122.4540 37.78251
## 3288  -122.4698 37.75298
## 3296  -122.4607 37.73104
## 3306  -122.4081 37.77962
## 3354  -122.4849 37.74663
## 3403  -122.4209 37.75792
## 3436  -122.4854 37.74788
## 3471  -122.4856 37.78116
## 3485  -122.4177 37.80345
## 3516  -122.4354 37.72201
## 3578  -122.4199 37.77606
## 3611  -122.4269 37.80229
## 3658  -122.4108 37.76534
## 3754  -122.3888 37.72872
## 3772  -122.3891 37.73264
## 4023  -122.4827 37.77762
## 4025  -122.4450 37.76430
## 4058  -122.4328 37.71546
## 4076  -122.4660 37.74828
## 4111  -122.5023 37.74333
## 4114  -122.4970 37.74468
## 4140  -122.4741 37.73601
## 4150  -122.4060 37.79853
## 4249  -122.4302 37.74726
## 4270  -122.4194 37.76710
## 4287  -122.4676 37.71048
## 4369  -122.4950 37.73327
## 4376  -122.3919 37.78868
## 4393  -122.3972 37.77570
## 4427  -122.4332 37.77869
## 4466  -122.4786 37.77394
## 4469  -122.4210 37.76470
## 4483  -122.4318 37.71761
## 4514  -122.4214 37.78139
## 4526  -122.4153 37.74047
## 4527  -122.4209 37.79368
## 4555  -122.4145 37.79167
## 4561  -122.3910 37.78938
## 4566  -122.4031 37.79058
## 4578  -122.3906 37.75501
## 4579  -122.4154 37.74537
## 4592  -122.3922 37.78607
## 4599  -122.4385 37.75980
## 4602  -122.4785 37.76295
## 4618  -122.4137 37.79099
## 4648  -122.4635 37.76461
## 4663  -122.4160 37.74563
## 4669  -122.3918 37.77899
## 4759  -122.4911 37.75318
## 4800  -122.4137 37.79099
## 4825  -122.4147 37.72528
## 4853  -122.4388 37.75978
## 4864  -122.3696 37.72900
## 4896  -122.4245 37.79493
## 4899  -122.3902 37.75236
## 4921  -122.5103 37.77369
## 4939  -122.3930 37.73299
## 4944  -122.3919 37.78868
## 4970  -122.4646 37.77383
## 4986  -122.3904 37.78139
## 5053  -122.4287 37.77484
## 5085  -122.4244 37.78587
## 5104  -122.4255 37.72744
## 5108  -122.4364 37.77781
## 5109  -122.4394 37.80326
## 5114  -122.4853 37.75287
## 5188  -122.4302 37.74806
## 5228  -122.3881 37.73174
## 5247  -122.4242 37.80551
## 5252  -122.4296 37.79988
## 5261  -122.4048 37.79879
## 5287  -122.4335 37.79946
## 5288  -122.4611 37.73481
## 5316  -122.3918 37.77899
## 5320  -122.4481 37.77945
## 5321  -122.4477 37.71766
## 5334  -122.3919 37.78868
## 5350  -122.4336 37.72535
## 5351  -122.4540 37.78251
## 5360  -122.4264 37.72624
## 5366  -122.3922 37.78607
## 5398  -122.3905 37.78759
## 5425  -122.3919 37.78868
## 5430  -122.4263 37.72116
## 5450  -122.4175 37.79865
## 5461  -122.3989 37.78577
## 5467  -122.4175 37.77341
## 5511  -122.3919 37.78868
## 5519  -122.4093 37.74024
## 5521  -122.4816 37.77294
## 5522  -122.4194 37.75596
## 5547  -122.3959 37.79054
## 5561  -122.4213 37.78884
## 5571  -122.4194 37.76710
## 5609  -122.4245 37.79920
## 5647  -122.4085 37.75911
## 5654  -122.4673 37.75439
## 5671  -122.4459 37.75811
## 5691  -122.4806 37.77724
## 5705  -122.4167 37.76091
## 5712  -122.3702 37.72838
## 5732  -122.4207 37.73596
## 5740  -122.4101 37.77808
## 5757  -122.4752 37.71513
## 5765  -122.3922 37.73203
## 5827  -122.4036 37.71724
## 5978  -122.4264 37.73053
## 5997  -122.3886 37.77159
## 6043  -122.5068 37.77785
## 6056  -122.4089 37.79286
## 6134  -122.4191 37.70898
## 6140  -122.4107 37.76483
## 6181  -122.3974 37.78575
## 6191  -122.4321 37.73585
## 6208  -122.4774 37.75583
## 6252  -122.4355 37.74743
## 6258  -122.3904 37.78271
## 6342  -122.4846 37.75914
## 6365  -122.4356 37.75016
## 6410  -122.3824 37.72887
## 6430  -122.4699 37.75813
## 6477  -122.4656 37.71671
## 6581  -122.4792 37.77343
## 6687  -122.4243 37.79352
## 6691  -122.4371 37.71173
## 6694  -122.4399 37.71837
## 6721  -122.4308 37.73647
## 6723  -122.4566 37.75242
## 6782  -122.4298 37.71104
## 6787  -122.4405 37.80504
## 6814  -122.4848 37.73940
## 6863  -122.4122 37.79293
## 6892  -122.4419 37.74582
## 6903  -122.4084 37.79154
## 6915  -122.4144 37.77516
## 6918  -122.3958 37.77816
## 6932  -122.4244 37.74805
## 6933  -122.4368 37.78249
## 6988  -122.5064 37.77555
## 7006  -122.4528 37.76708
## 7009  -122.4947 37.76067
## 7032  -122.4050 37.77762
## 7073  -122.4194 37.75596
## 7142  -122.4953 37.73096
## 7149  -122.4104 37.75357
## 7175  -122.4370 37.76934
## 7206  -122.4936 37.77938
## 7303  -122.4327 37.77299
## 7320  -122.4194 37.76710
## 7337  -122.4752 37.71551
## 7395  -122.4263 37.76856
## 7404  -122.3939 37.75915
## 7447  -122.4113 37.72850
## 7462  -122.4127 37.76517
## 7493  -122.4194 37.75596
## 7533  -122.4515 37.75418
## 7555  -122.4071 37.77553
## 7564  -122.4389 37.73129
## 7600  -122.4097 37.77613
## 7649  -122.4290 37.78872
## 7651  -122.3958 37.77816
## 7704  -122.4337 37.76214
## 7735  -122.3904 37.78271
## 7832  -122.4853 37.75355
## 7833  -122.3983 37.75753
## 7848  -122.4490 37.76795
## 7936  -122.4087 37.77310
## 7981  -122.4489 37.71967
## 7991  -122.4325 37.77805
## 8047  -122.4019 37.78005
## 8089  -122.4194 37.75596
## 8127  -122.4685 37.76188
## 8161  -122.4307 37.74163
## 8275  -122.4126 37.77870
## 8322  -122.3924 37.73700
## 8344  -122.4126 37.77870
## 8378  -122.4900 37.78324
## 8401  -122.4642 37.71638
## 8410  -122.4279 37.77408
## 8419  -122.4016 37.75145
## 8422  -122.4836 37.71006
## 8438  -122.4668 37.73426
## 8517  -122.3919 37.78868
## 8537  -122.3980 37.71952
## 8590  -122.4291 37.73289
## 8612  -122.4194 37.75596
## 8652  -122.4222 37.72570
## 8694  -122.3919 37.78868
## 8703  -122.3958 37.77339
## 8713  -122.4301 37.77286
## 8730  -122.4126 37.77870
## 8818  -122.3897 37.78278
## 8823  -122.4361 37.78937
## 8842  -122.4257 37.79707
## 8844  -122.4259 37.78673
## 8888  -122.4318 37.79685
## 8891  -122.4277 37.77402
## 8899  -122.3896 37.78063
## 8907  -122.4728 37.75198
## 8932  -122.4102 37.74514
## 8962  -122.3942 37.77181
## 8983  -122.4206 37.80236
## 9049  -122.4651 37.75119
## 9054  -122.4371 37.73361
## 9179  -122.4262 37.72086
## 9216  -122.4160 37.73419
## 9219  -122.4550 37.74285
## 9356  -122.4522 37.73356
## 9376  -122.3982 37.75309
## 9384  -122.4152 37.73405
## 9469  -122.4174 37.80442
## 9479  -122.4251 37.72966
## 9490  -122.4835 37.74208
## 9569  -122.4723 37.73269
## 9639  -122.4883 37.77411
## 9640  -122.4509 37.73126
## 9647  -122.3909 37.75967
## 9675  -122.3891 37.76366
## 9691  -122.4381 37.78210
## 9764  -122.4875 37.76255
## 9776  -122.4422 37.72376
## 9806  -122.3962 37.77393
## 9830  -122.4642 37.71673
## 9906  -122.4433 37.76135
## 9959  -122.3934 37.77764
## 9966  -122.4639 37.72971
## 9985  -122.4411 37.75920
## 10010 -122.4279 37.76639
## 10046 -122.4755 37.71495
## 10063 -122.3879 37.73919
## 10066 -122.4429 37.78894
## 10100 -122.3890 37.73260
## 10157 -122.4377 37.79704
## 10164 -122.4180 37.74667
## 10170 -122.4453 37.79137
## 10229 -122.4137 37.79099
## 10239 -122.4451 37.78484
## 10257 -122.4102 37.74590
## 10263 -122.4353 37.75147
## 10272 -122.4974 37.78118
## 10273 -122.4401 37.75541
## 10307 -122.4167 37.79138
## 10338 -122.4343 37.73743
## 10343 -122.4108 37.77597
## 10345 -122.5037 37.74846
## 10351 -122.3910 37.78938
## 10353 -122.4865 37.76343
## 10410 -122.3940 37.75300
## 10415 -122.4287 37.73629
## 10419 -122.4344 37.76357
## 10499 -122.4531 37.73537
## 10530 -122.4150 37.73959
## 10578 -122.4594 37.73954
## 10583 -122.4760 37.71558
## 10600 -122.4832 37.78375
## 10607 -122.4911 37.73160
## 10632 -122.4347 37.78713
## 10646 -122.4518 37.76347
## 10661 -122.4066 37.73036
## 10781 -122.4124 37.80434
## 10782 -122.4554 37.71773
## 10800 -122.4152 37.73705
## 10857 -122.4178 37.76399
## 10887 -122.4260 37.74559
## 10914 -122.4349 37.74579
## 10957 -122.4173 37.74421
## 10967 -122.4587 37.73858
## 10985 -122.4325 37.77540
## 11005 -122.4558 37.73265
## 11023 -122.4173 37.73753
## 11035 -122.4386 37.76810
## 11094 -122.4275 37.72671
## 11158 -122.3937 37.78702
## 11160 -122.4236 37.77800
## 11170 -122.5054 37.76363
## 11171 -122.4001 37.79803
## 11180 -122.4757 37.71556
## 11181 -122.4359 37.75531
## 11202 -122.4133 37.78070
## 11258 -122.4437 37.80053
## 11265 -122.4865 37.76465
## 11284 -122.4117 37.75704
## 11286 -122.4194 37.75596
## 11305 -122.4136 37.79298
## 11385 -122.4424 37.74562
## 11451 -122.4370 37.73761
## 11498 -122.4438 37.71852
## 11554 -122.4990 37.77913
## 11564 -122.4194 37.75596
## 11576 -122.3910 37.78938
## 11598 -122.3824 37.72852
## 11603 -122.4748 37.71490
## 11622 -122.3988 37.71372
## 11636 -122.4678 37.72225
## 11690 -122.4415 37.78044
## 11728 -122.4793 37.75869
## 11729 -122.4165 37.73373
## 11808 -122.4457 37.75769
## 11830 -122.4965 37.77600
## 11844 -122.4468 37.73736
## 11873 -122.4474 37.73477
## 11876 -122.4126 37.79035
## 11882 -122.3970 37.75975
## 11918 -122.4438 37.76155
## 11927 -122.4333 37.73432
## 11968 -122.3936 37.78370
## 11969 -122.3907 37.78117
## 11971 -122.4742 37.75636
## 11972 -122.4542 37.77443
## 12052 -122.4385 37.75442
## 12232 -122.5060 37.73665
## 12285 -122.4786 37.77955
## 12350 -122.4672 37.73270
## 12383 -122.4093 37.77680
## 12428 -122.4274 37.77132
## 12467 -122.3971 37.78091
## 12533 -122.4024 37.77546
## 12559 -122.4130 37.80397
## 12569 -122.4114 37.72936
## 12622 -122.5059 37.73682
## 12668 -122.4131 37.79777
## 12672 -122.4154 37.73877
## 12675 -122.4357 37.73740
## 12730 -122.4011 37.73263
## 12770 -122.4221 37.72534
## 12783 -122.3910 37.78938
## 12865 -122.3919 37.78868
## 12866 -122.4891 37.75507
## 13000 -122.3968 37.77945
## 13015 -122.4341 37.76397
## 13054 -122.4165 37.77528
## 13177 -122.4043 37.78622
## 13251 -122.4294 37.77067
## 13252 -122.3933 37.77891
## 13268 -122.4316 37.76576
## 13293 -122.4936 37.77938
## 13324 -122.4116 37.75238
## 13382 -122.4513 37.76809
## 13391 -122.4248 37.80057
## 13414 -122.4433 37.78314
## 13422 -122.4141 37.80393
## 13470 -122.4439 37.71527
## 13472 -122.4815 37.74502
## 13478 -122.4214 37.78139
## 13507 -122.4298 37.73805
## 13529 -122.4242 37.80551
## 13550 -122.4439 37.80384
## 13581 -122.3956 37.73468
## 13593 -122.4286 37.74623
## 13595 -122.4617 37.78577
## 13662 -122.3918 37.77899
## 13675 -122.3696 37.72900
## 13701 -122.3962 37.77393
## 13722 -122.4025 37.73527
## 13777 -122.4750 37.71499
## 13821 -122.4032 37.78824
## 13895 -122.4485 37.73998
## 13903 -122.3919 37.78868
## 13926 -122.4141 37.75915
## 13936 -122.4296 37.73278
## 14044 -122.4484 37.76533
## 14072 -122.4546 37.74500
## 14087 -122.4642 37.75929
## 14239 -122.4665 37.71392
## 14241 -122.4456 37.80221
## 14256 -122.4990 37.77356
## 14287 -122.4234 37.77253
## 14338 -122.4475 37.77957
## 14345 -122.4192 37.75400
## 14353 -122.3696 37.72900
## 14367 -122.4253 37.79247
## 14413 -122.4397 37.76158
## 14496 -122.4276 37.79552
## 14517 -122.4617 37.75037
## 14679 -122.4565 37.72132
## 14690 -122.4365 37.71167
## 14699 -122.4375 37.76305
## 14721 -122.3919 37.78868
## 14728 -122.3934 37.77764
## 14769 -122.4258 37.78636
## 14773 -122.4326 37.75640
## 14777 -122.3936 37.78370
## 14803 -122.4121 37.77255
## 14870 -122.4481 37.73958
## 14965 -122.4093 37.74690
## 14974 -122.5070 37.78027
## 15069 -122.4340 37.74511
## 15076 -122.4733 37.77618
## 15143 -122.4234 37.77253
## 15159 -122.4937 37.78119
## 15167 -122.4129 37.78791
## 15179 -122.4668 37.74374
## 15186 -122.4476 37.77794
## 15187 -122.3991 37.78318
## 15197 -122.3943 37.77898
## 15208 -122.4557 37.71642
## 15217 -122.4678 37.72282
## 15260 -122.4416 37.71429
## 15267 -122.4507 37.74389
## 15319 -122.4345 37.78565
## 15327 -122.4957 37.77863
## 15328 -122.4564 37.70934
## 15349 -122.4151 37.73756
## 15357 -122.3761 37.73193
## 15426 -122.4009 37.76461
## 15474 -122.4009 37.76461
## 15485 -122.5046 37.75155
## 15486 -122.4718 37.76104
## 15498 -122.4174 37.74938
## 15525 -122.4124 37.77093
## 15551 -122.3943 37.78613
## 15567 -122.4055 37.80424
## 15580 -122.3914 37.78451
## 15587 -122.5046 37.74607
## 15589 -122.4339 37.72119
## 15604 -122.4248 37.75359
## 15615 -122.4098 37.71496
## 15621 -122.4667 37.73841
## 15666 -122.4379 37.75930
## 15669 -122.4194 37.75596
## 15713 -122.3702 37.72838
## 15717 -122.4382 37.75646
## 15790 -122.4021 37.76545
## 15887 -122.4311 37.76352
## 15918 -122.3942 37.77181
## 15933 -122.4178 37.76399
## 15945 -122.4390 37.79617
## 16033 -122.4014 37.78629
## 16076 -122.4149 37.79362
## 16105 -122.4489 37.76179
## 16173 -122.4135 37.72849
## 16199 -122.4014 37.77974
## 16204 -122.4327 37.76614
## 16210 -122.4210 37.75876
## 16212 -122.4126 37.77870
## 16215 -122.3958 37.77705
## 16229 -122.4400 37.71831
## 16310 -122.3988 37.72731
## 16315 -122.4254 37.80339
## 16389 -122.4800 37.75339
## 16440 -122.4970 37.74495
## 16450 -122.4194 37.75596
## 16458 -122.4624 37.77409
## 16474 -122.4542 37.77443
## 16480 -122.4014 37.77974
## 16500 -122.4126 37.77870
## 16507 -122.4300 37.78804
## 16517 -122.4240 37.76628
## 16531 -122.4266 37.73311
## 16549 -122.5038 37.74097
## 16640 -122.4384 37.77517
## 16643 -122.4517 37.74504
## 16712 -122.4194 37.75596
## 16755 -122.3895 37.73005
## 16778 -122.4708 37.75810
## 16796 -122.4035 37.72002
## 16800 -122.3696 37.72900
## 16802 -122.4042 37.71115
## 16921 -122.3696 37.72900
## 16934 -122.4054 37.77794
## 17006 -122.4459 37.77465
## 17051 -122.4560 37.75650
## 17073 -122.4386 37.73261
## 17084 -122.3942 37.77181
## 17091 -122.4750 37.71499
## 17156 -122.4286 37.71810
## 17168 -122.4766 37.76075
## 17199 -122.3910 37.78938
## 17200 -122.4151 37.78524
## 17223 -122.4011 37.77047
## 17249 -122.4159 37.74769
## 17252 -122.4908 37.73143
## 17265 -122.4043 37.80443
## 17304 -122.4014 37.77974
## 17345 -122.4245 37.79493
## 17358 -122.4704 37.72549
## 17384 -122.3904 37.78271
## 17404 -122.4080 37.77576
## 17410 -122.3904 37.75368
## 17465 -122.5033 37.74841
## 17491 -122.4015 37.73255
## 17510 -122.4359 37.75497
## 17517 -122.4519 37.76025
## 17542 -122.4124 37.77093
## 17562 -122.4514 37.70881
## 17613 -122.3893 37.73988
## 17618 -122.4249 37.75375
## 17698 -122.4262 37.79417
## 17715 -122.4140 37.75590
## 17754 -122.4360 37.80400
## 17767 -122.4916 37.78122
## 17857 -122.4065 37.77812
## 17867 -122.4777 37.76001
## 17876 -122.5107 37.77399
## 17881 -122.4226 37.73575
## 17887 -122.4188 37.77159
## 17912 -122.3886 37.77159
## 17929 -122.4289 37.72161
## 17933 -122.4452 37.71900
## 17936 -122.4458 37.72772
## 17986 -122.4429 37.73288
## 17988 -122.4473 37.73331
## 18002 -122.4422 37.77803
## 18018 -122.4772 37.74409
## 18050 -122.4198 37.76161
## 18053 -122.4038 37.76405
## 18059 -122.4221 37.74746
## 18071 -122.4271 37.73809
## 18107 -122.4388 37.73398
## 18164 -122.4160 37.73504
## 18227 -122.4805 37.75531
## 18231 -122.4671 37.76370
## 18323 -122.4229 37.80586
## 18363 -122.3727 37.73057
## 18386 -122.4767 37.77714
## 18433 -122.4408 37.78769
## 18439 -122.4345 37.73587
## 18444 -122.4716 37.78657
## 18453 -122.3928 37.78205
## 18470 -122.4682 37.75544
## 18504 -122.5075 37.74854
## 18511 -122.3959 37.79054
## 18553 -122.3942 37.77181
## 18584 -122.4140 37.75590
## 18619 -122.4315 37.73498
## 18627 -122.4716 37.75516
## 18642 -122.4780 37.78258
## 18653 -122.4258 37.72884
## 18718 -122.4391 37.75547
## 18733 -122.4294 37.72902
## 18742 -122.4194 37.75596
## 18792 -122.4181 37.78824
## 18795 -122.3938 37.78638
## 18801 -122.4288 37.72584
## 18871 -122.4234 37.77253
## 18886 -122.4305 37.74402
## 18913 -122.4280 37.74869
## 18960 -122.4387 37.71054
## 18968 -122.4041 37.80347
## 18971 -122.4058 37.80159
## 18990 -122.4398 37.75879
## 19015 -122.4613 37.76457
## 19065 -122.4121 37.77187
## 19082 -122.4684 37.76042
## 19105 -122.4836 37.71006
## 19143 -122.4642 37.78458
## 19153 -122.4520 37.73468
## 19197 -122.4102 37.75679
## 19202 -122.4412 37.71904
## 19227 -122.3896 37.78063
## 19281 -122.4309 37.72308
## 19325 -122.4287 37.77642
## 19375 -122.4152 37.73366
## 19532 -122.4596 37.78655
## 19535 -122.5019 37.77481
## 19570 -122.4623 37.73544
## 19626 -122.4348 37.74426
## 19641 -122.3919 37.78868
## 19718 -122.4194 37.75596
## 19765 -122.4194 37.75596
## 19782 -122.3911 37.73912
## 19802 -122.4152 37.79845
## 19811 -122.3951 37.72508
## 19832 -122.4798 37.78437
## 19838 -122.4301 37.74120
## 19845 -122.4435 37.75964
## 19880 -122.4243 37.79352
## 19897 -122.3934 37.77764
## 19932 -122.4590 37.71972
## 19933 -122.4261 37.73585
## 19954 -122.4934 37.76073
## 19991 -122.4242 37.79788
## 20009 -122.4528 37.71555
## 20026 -122.3909 37.75967
## 20061 -122.4382 37.80320
## 20095 -122.3919 37.78868
## 20102 -122.4936 37.74282
## 20221 -122.4251 37.75560
## 20235 -122.4709 37.75704
## 20293 -122.4937 37.77999
## 20310 -122.4252 37.74835
## 20320 -122.4355 37.71758
## 20321 -122.4739 37.74485
## 20347 -122.5013 37.76014
## 20356 -122.4517 37.73356
## 20426 -122.4014 37.77974
## 20484 -122.4345 37.80458
## 20562 -122.3974 37.78575
## 20605 -122.4729 37.77471
## 20611 -122.3922 37.78607
## 20655 -122.4621 37.77642
## 20687 -122.4548 37.78306
## 20743 -122.4469 37.71901
## 20754 -122.4493 37.77461
## 20768 -122.3962 37.77393
## 20778 -122.4743 37.76473
## 20779 -122.4024 37.76117
## 20804 -122.5022 37.74887
## 20871 -122.3922 37.73197
## 20877 -122.4180 37.76654
## 20928 -122.4731 37.77999
## 20971 -122.4816 37.77795
## 20979 -122.4582 37.73860
## 21007 -122.3696 37.72900
## 21053 -122.4194 37.75596
## 21064 -122.4263 37.76856
## 21079 -122.4425 37.80069
## 21093 -122.4336 37.74361
## 21105 -122.4245 37.79920
## 21124 -122.4399 37.79230
## 21152 -122.4642 37.72043
## 21174 -122.4227 37.79931
## 21201 -122.4340 37.77121
## 21231 -122.4272 37.72286
## 21255 -122.4094 37.75902
## 21374 -122.4208 37.73601
## 21424 -122.3702 37.72838
## 21455 -122.4195 37.80612
## 21456 -122.3932 37.78006
## 21458 -122.4980 37.75912
## 21479 -122.4331 37.79233
## 21486 -122.3922 37.78607
## 21496 -122.4739 37.78414
## 21499 -122.4031 37.79058
## 21508 -122.4606 37.72592
## 21522 -122.4188 37.77159
## 21599 -122.4630 37.74288
## 21626 -122.4271 37.78448
## 21632 -122.4368 37.74525
## 21644 -122.4176 37.76271
## 21652 -122.3901 37.75370
## 21678 -122.3702 37.72838
## 21694 -122.4414 37.79054
## 21711 -122.4126 37.77870
## 21712 -122.3904 37.78139
## 21790 -122.4446 37.79053
## 21791 -122.5024 37.73799
## 21834 -122.4218 37.79478
## 21889 -122.4157 37.79792
## 21932 -122.4520 37.71732
## 21957 -122.4282 37.79869
## 21972 -122.4552 37.74086
## 21985 -122.4549 37.73345
## 22007 -122.4195 37.74526
## 22024 -122.5006 37.75075
## 22064 -122.4569 37.75564
## 22085 -122.4497 37.72884
## 22111 -122.3933 37.77891
## 22167 -122.5041 37.73937
## 22254 -122.3992 37.79825
## 22255 -122.4107 37.77166
## 22273 -122.4448 37.71416
## 22324 -122.3986 37.71395
## 22340 -122.4376 37.75052
## 22499 -122.4438 37.72965
## 22531 -122.4141 37.77608
## 22533 -122.4409 37.74887
## 22549 -122.4431 37.71275
## 22554 -122.4009 37.76461
## 22600 -122.4362 37.74701
## 22607 -122.4753 37.71516
## 22673 -122.4131 37.74787
## 22739 -122.3919 37.78868
## 22801 -122.4709 37.75504
## 22802 -122.4014 37.77974
## 22820 -122.4887 37.77819
## 22870 -122.4367 37.76372
## 22871 -122.4507 37.72158
## 22896 -122.4342 37.74955
## 22924 -122.4460 37.71426
## 22931 -122.4210 37.74497
## 22951 -122.4948 37.77410
## 23000 -122.4254 37.76334
## 23036 -122.4163 37.80401
## 23045 -122.3996 37.77438
## 23056 -122.4753 37.71487
## 23074 -122.3696 37.72900
## 23117 -122.4223 37.79267
## 23135 -122.3942 37.77181
## 23154 -122.4648 37.75510
## 23172 -122.4415 37.78044
## 23173 -122.4835 37.78257
## 23251 -122.4119 37.71432
## 23267 -122.4263 37.71022
## 23308 -122.4769 37.76468
## 23356 -122.4613 37.71719
## 23400 -122.4689 37.77570
## 23433 -122.4753 37.71553
## 23467 -122.3928 37.78205
## 23478 -122.4855 37.75523
## 23486 -122.4276 37.72233
## 23503 -122.4108 37.76534
## 23519 -122.4454 37.76390
## 23535 -122.4639 37.76393
## 23540 -122.4599 37.72230
## 23542 -122.4860 37.78680
## 23598 -122.3910 37.78938
## 23601 -122.3922 37.78607
## 23721 -122.4452 37.71600
## 23782 -122.4012 37.75378
## 23790 -122.5038 37.75594
## 23828 -122.3957 37.73636
## 23881 -122.3906 37.75501
## 23895 -122.4517 37.72883
## 23896 -122.4011 37.77047
## 23913 -122.4734 37.75749
## 23918 -122.4676 37.71048
## 23919 -122.4936 37.77938
## 23991 -122.4418 37.73477
## 24130 -122.4544 37.76267
## 24185 -122.4891 37.77817
## 24223 -122.4194 37.75596
## 24282 -122.4934 37.73923
## 24314 -122.4577 37.74434
## 24326 -122.3880 37.73296
## 24374 -122.3696 37.72900
## 24393 -122.4372 37.75934
## 24401 -122.4517 37.73329
## 24402 -122.4711 37.74663
## 24504 -122.4536 37.70915
## 24546 -122.4544 37.78545
## 24554 -122.4395 37.75941
## 24574 -122.3997 37.75216
## 24607 -122.4836 37.77285
## 24639 -122.4824 37.75738
## 24659 -122.3920 37.72865
## 24674 -122.4760 37.78306
## 24700 -122.4138 37.73519
## 24716 -122.4922 37.77477
## 24724 -122.4374 37.75893
## 24742 -122.4704 37.78430
## 24781 -122.4137 37.79099
## 24794 -122.4016 37.75125
## 24816 -122.4009 37.72799
## 24840 -122.4704 37.74936
## 24872 -122.4236 37.77800
## 24896 -122.4441 37.73437
## 24922 -122.4446 37.77193
## 24940 -122.4192 37.72509
## 24949 -122.4356 37.72596
## 24966 -122.3919 37.78868
## 24983 -122.4458 37.73512
## 24986 -122.4668 37.77558

str(sfhomes15) # the data frame

str(sfhomes15_sp) # the data frame

str(sfhomes15) # the data frame

## 'data.frame':    835 obs. of  19 variables:
##  $ FiscalYear      : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ SalesDate       : Factor w/ 5389 levels "1920-01-30","1920-01-31",..: 5304 5298 5387 5307 5270 5256 5336 5216 5283 5323 ...
##  $ Address         : Factor w/ 24696 levels "0000 0000                       0000",..: 19602 9677 10982 3007 12212 14195 4416 16847 6820 21871 ...
##  $ YearBuilt       : int  1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
##  $ NumBedrooms     : int  2 2 2 3 3 3 3 2 0 2 ...
##  $ NumBathrooms    : int  2 2 2 3 3 3 3 2 0 2 ...
##  $ NumRooms        : int  5 5 5 6 7 7 7 5 0 0 ...
##  $ NumStories      : int  0 1 1 0 1 14 2 1 0 1 ...
##  $ NumUnits        : int  1 1 1 1 1 1 1 1 0 1 ...
##  $ AreaSquareFeet  : int  1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
##  $ ImprovementValue: int  432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
##  $ LandValue       : int  432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
##  $ Neighborhood    : Factor w/ 40 levels "","Bayview Hunters Point",..: 39 29 20 14 19 21 8 27 14 7 ...
##  $ Location        : Factor w/ 20416 levels "(37.7081695281752, -122.449582481739)",..: 5483 11922 14534 1141 11757 18906 6555 18389 1091 18367 ...
##  $ SupeDistrict    : int  7 10 6 7 9 3 8 2 7 6 ...
##  $ totvalue        : int  865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
##  $ SalesYear       : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ lat             : num  37.7 37.8 37.8 37.7 37.8 ...
##  $ lon             : num  -122 -122 -122 -122 -122 ...

str(sfhomes15_sp) # the SPDF

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  835 obs. of  17 variables:
##   .. ..$ FiscalYear      : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   .. ..$ SalesDate       : Factor w/ 5389 levels "1920-01-30","1920-01-31",..: 5304 5298 5387 5307 5270 5256 5336 5216 5283 5323 ...
##   .. ..$ Address         : Factor w/ 24696 levels "0000 0000                       0000",..: 19602 9677 10982 3007 12212 14195 4416 16847 6820 21871 ...
##   .. ..$ YearBuilt       : int [1:835] 1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
##   .. ..$ NumBedrooms     : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
##   .. ..$ NumBathrooms    : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
##   .. ..$ NumRooms        : int [1:835] 5 5 5 6 7 7 7 5 0 0 ...
##   .. ..$ NumStories      : int [1:835] 0 1 1 0 1 14 2 1 0 1 ...
##   .. ..$ NumUnits        : int [1:835] 1 1 1 1 1 1 1 1 0 1 ...
##   .. ..$ AreaSquareFeet  : int [1:835] 1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
##   .. ..$ ImprovementValue: int [1:835] 432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
##   .. ..$ LandValue       : int [1:835] 432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
##   .. ..$ Neighborhood    : Factor w/ 40 levels "","Bayview Hunters Point",..: 39 29 20 14 19 21 8 27 14 7 ...
##   .. ..$ Location        : Factor w/ 20416 levels "(37.7081695281752, -122.449582481739)",..: 5483 11922 14534 1141 11757 18906 6555 18389 1091 18367 ...
##   .. ..$ SupeDistrict    : int [1:835] 7 10 6 7 9 3 8 2 7 6 ...
##   .. ..$ totvalue        : int [1:835] 865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
##   .. ..$ SalesYear       : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   ..@ coords.nrs : int [1:2] 19 18
##   ..@ coords     : num [1:835, 1:2] -122 -122 -122 -122 -122 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:835] "24" "35" "76" "79" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] -122.5 37.7 -122.4 37.8
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

You can see from str(sfhomes) that a SPDF object is a collection of slots or components. The key ones are:

• @data data frame of attributes that describe each location
• @coords the coordinates for each geometric object - here points
• @bbox the min and max bounding coordinates
• @proj4string the coordinate reference system defintion as a string

Explore the object in the Environment window

Review the output of each of these:

summary(sfhomes15_sp)
head(sfhomes15_sp@data)
class(sfhomes15_sp@data)

sfhomes15_sp@bbox
bbox(sfhomes15_sp)

head(sfhomes15_sp@coords)
head(sfhomes15_sp$lat)
head(sfhomes15_sp$lon)

sfhomes15_sp@proj4string
proj4string(sfhomes15_sp)

Look at sfhomes15_sp in the environment window

Are all the columns that were present in sfhomes15 also in sfhomes15_sp?

Is there a slot in sfhomes15_sp without data?

proj4string(sfhomes15_sp) # get a CRS object

## [1] NA

plot(sfhomes15_sp)  # using sp::plot

We created great maps of sfhomes point data with ggplot and ggmap.

Then we created a simple map of the SPDF with plot.

We aren't seeing the value of sp objects just yet.

There's an R package for that!

rgdal is an R port of the powerful and widely used GDAL library.

It is the most commonly used R library for importing and exporting spatial data.

• OGR: for vector data: readOGR() and writeOGR()

• GDAL for raster data: readGDAL() and writeGDAL()

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.4.3

## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-5

# See what file types are supported by rgdal drivers
# ogrDrivers()$name

gdal.org

?readOGR

For more info on working with rgdal to load different types of spatial data in R see this excellent tutorial by Zev Ross.

Take a look at the file(s)

dir("data", pattern="sf_boundary")

## [1] "sf_boundary.dbf" "sf_boundary.prj" "sf_boundary.shp" "sf_boundary.shx"

sfboundary <- readOGR(dsn="data",layer="sf_boundary")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_boundary"
## with 1 features
## It has 3 fields

# or
# sfboundary <- readOGR("data","sf_boundary")
# but not
#sfboundary <- readOGR(dsn="data/",layer="sf_boundary")

What type of sp object is sfboundary? How many features are in the object?

class(sfboundary)
str(sfboundary) 
head(sfboundary@data)  

Explore the object in the Envi window

How?

plot(sfboundary)

Is it a geographic or projected CRS?

proj4string(sfboundary)

plot(sfboundary)
points(sfhomes15_sp, col="red")

Where are the points? What's wrong?

plot(sfboundary)
points(sfhomes15_sp, col="red")

Compare the CRSs, are they the same?

proj4string(sfboundary)
proj4string(sfhomes15_sp)
proj4string(sfboundary) == proj4string(sfhomes15_sp)

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] NA

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] NA

sfboundary@bbox

##         min       max
## x  542696.6  556659.9
## y 4173563.7 4185088.6

sfhomes15_sp@bbox

##            min        max
## lon -122.51072 -122.36964
## lat   37.70854   37.80612

The #1 reason…

All sp objects should have a defined CRS

If not, one needs to be assigned to the object

• This is called defining a projection

• This doesn't change the coordinates

All sp objects should have the same CRS.

• When they don't, they need to be transformed to a common CRS.

• This is called a projection transformation,
  • or projecting or reprojection.

• Projection transformation returns a new spatial object with the transformed coordinates

GIS software - or software for working with geospatial data - will have a database of definitions for thousands of Earth referenced coordinate systems

and methods for using those CRS definitions to define or transform a CRS.

So, do we need to define the CRS of either

• sfboundary or

• sfhomes15_sp?

If yes, which one? How can you tell?

We need to know

• what the appropriate CRS is
• how to define it as a CRS object
• how to assign it to an sp object

bbox(sfhomes15_sp)

##            min        max
## lon -122.51072 -122.36964
## lat   37.70854   37.80612

Use the sp function CRS() to define a coordinate reference system

• requires as input a string that defines the CRS parameters

Then, use the sp function proj4string() to assign it

proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
                               +ellps=WGS84 +datum=WGS84 +no_defs")  

The WGS84 geographic coordinate referenece system can be defined by:

1. It's paramater string

CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

1. A code that can be used to look up the parameters

CRS("+init=epsg:4326")

EPSG refers to the group that defined these

• European Petroleum Survey Group

EPSG codes are more commonly used. Can you guess why?

# use an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326") 

# or enter the parameter string
# proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
#                               +ellps=WGS84 +datum=WGS84 +no_defs")  

Proj4 is the standard library for defining and transforming map projections and coordinate reference systems.

Here is an example file of proj4 CRS definitions

Once you set the CRS you can check it with the same function.

proj4string(sfhomes15_sp)

## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

proj4string() can get or set the CRS

CRS() is used to define the CRS of an sp object

The following are equivalent

# Define the CRS with the parameter string
proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
                              +ellps=WGS84 +datum=WGS84 +no_defs")  

# Define the CRS with an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326") 

epsg codes are used as short-hand for CRSs definitions

Are they both defined? Equal?

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] FALSE

Use sp function spTransform

Requires as input:

• a sp object to transform with a defined CRS

• a target CRS

Outputs a new spatial object with coordinate data in the target CRS

• It does not change the input data

Let's transform the CRS of sfboundary to be the same as that of sfhomes15_sp.

Why not the other way around?

sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))

Do they match now?

How will we know?

proj4string(sfhomes15_sp) == proj4string(sfboundary_lonlat)

## [1] FALSE

plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")

Woo-hoo!

We can transform sp objects to the same CRS for mapping and spatial analysis (we'll get there!)

Use writeOGR to save sfboundary_lonlat to a new shapefile

See ?writeOGR for help

# write transformed data to a new shapefile 
writeOGR(sfboundary_lonlat, 
          dsn = "data", 
          layer = "sfbounary_lonlat", 
          driver="ESRI Shapefile")

# is it there?
dir("data")

We want all data in the same CRS

Which one is best?

Geographic CRSs

• 4326 Geographic, WGS84 (default for lon/lat)

• 4269 Geographic, NAD83 (USA Fed agencies like Census)

Projected CRSs

• 5070 USA Contiguous Albers Equal Area Conic

• 3310 CA ALbers Equal Area

• 26910 UTM Zone 10, NAD83 (Northern Cal)

• 3857 Web Mercator (web maps)

See http://spatialreference.org/

Use this site to find EPSG codes and proj4 CRS strings

1. Use coordinates to convert the landmarks data frame to a SpatialPointsDataFrame
2. Use proj4string to define it's CRS as 3857 Web Mercator
3. Use spTransform to transform the landmarks data to geographic coords (WGS84)
4. Add the landmarks to the last map we created.

# Convert the DF to a SPDF
coordinates(landmarks) <-c("X","Y")
# Define the CRS with an EPSG code for Web mercator
proj4string(landmarks) <- CRS("+init=epsg:3857") 
# Transform the CRS to WGS84
landmarks_lonlat <- spTransform(landmarks, CRS("+init=epsg:4326"))
# map it
plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")

1. Use readOGR to load the shapefile "data/sf_highways.shp"
2. Check its data type with class()
3. Check its CRS with proj4string()
4. Use spTransform if need to project it to WGS84
5. Recreate the last map and add the highways as black lines

highways <- readOGR(dsn="data", layer="sf_highways")
class(highways)
proj4string(highways)
highways_lonlat <- spTransform(highways, CRS("+init=epsg:4326"))

plot(sfboundary_lonlat)
lines(highways_lonlat, col="black")
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_highways"
## with 246 features
## It has 5 fields

## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"

## [1] "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

# QUESTIONS?

base::plot, ggplot, ggmap for geospatial data in data frames

• great for creating maps given these types of data

sp::plot for sp objects

• meh, but great for a quick look at spatial data

Use it to create quick maps

Can be used to create great maps

• BUT complex, non-intuitive syntax = long ugly code

See examples in sp Gallery: Plotting maps with sp

# map of the sfhomes data by totalvaue
spplot(sfhomes15_sp,"totvalue")

tmap stands for thematic map

Great maps with less code than the alternatives

Syntax should be familar to ggplot2 users, but simpler

Relatively easy to create interactive maps

• ?tmap
• vignette("tmap-nutshell")
• vignette("tmap-modes")
• tmap in a Nutshell

Load the library

library(tmap)

## Warning: package 'tmap' was built under R version 3.4.3

qtm(sfhomes15_sp)

tmap_mode("view")

## tmap mode set to interactive viewing

qtm(sfhomes15_sp)

tmap_mode("plot")

## tmap mode set to plotting

ttm()       # Toggle Tmap mode between "interactive" and "plot"

last_map()  # display the last map 

tmap_mode("plot")

## tmap mode set to plotting

tmap's flexibility comes in how it intitively allows you to layer spatial data and style the layers by data attributes

Use tm_shape(<sp_object>) to specifiy a geospatial data layer

Add + tm_<element>(...) to style the layer by data values

…and other options for creating a publication ready map

?tmap_shape

?tmap_element

• ?tm_polygons (tm_fill, tm_borders)

• ?tm_symbols (tm_dots, etc…)

• tmap in a Nutshell

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black")

tm_shape(sfhomes15_sp) + 
  tm_dots(col="red", size=.25)

tm_shape(highways_lonlat) + 
  tm_lines(col="black")

## Mapping values

tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25)    # column names must be quoted

tm_dots are a type of tm_symbols()

See ?tm_symbols for other types of point symbols.

tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25)  # columnn names **must** be quoted!

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25) 

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25, title = "San Francisco Property Values (2015)") + 
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

Add the landmarks to the previous tmap

Feel free to change the colors and title etc.

What does the landmarks layer tell you about the way tmap handles CRSs?

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25, 
          title = "San Francisco Property Values (2015)") + 
tm_shape(landmarks) +
  tm_markers(col="black", size=.4, text="name",text.just = "left") +
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

Make the previous tmap interactive

ttm()       # Toggle Tmap mode between "interactive" and "plot"
last_map()  # display the last map 

## tmap mode set to interactive viewing

popup.vars = c("SalesYear","totvalue","NumBedrooms", "NumBathrooms","AreaSquareFeet")

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25, 
          title = "San Francisco Property Values (2015)",
          popup.vars=c("SalesYear","totvalue","NumBedrooms",
                       "NumBathrooms","AreaSquareFeet")) + 
tm_shape(landmarks) +
  tm_markers(col="black", size=.4, text="name",text.just = "left") +
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

Save it as an R object

map1 <- last_map()
map1 # then display it

Try this and then take a look at your files

save_tmap(map1, "sf_properties.png", height=6) # Static image file
save_tmap(map1, "sf_properties.html") # interactive web map

## Map saved to /Users/pattyf/Documents/Dlab/workshops/2018/r-geospatial-workshop/sf_properties.png

## Resolution: 2318.548 by 1800 pixels

## Size: 7.728492 by 6 inches (300 dpi)

## Interactive map saved to /Users/pattyf/Documents/Dlab/workshops/2018/r-geospatial-workshop/sf_properties.html

Many ways to do this.

You can share you map online by publishing it to RPubs.

• You need to have an RPubs account to make that work.

1. Enter the name of your tmap object (map1) or your tmap code in the console

2. In the Viewer window, click on the Publish icon.

tmap_mode("interactive")
map1 #  in console

http://rpubs.com/pfrontiera/358110

• Geospatial Data in R
• CSV > Data Frame > ggplot/ggmap
• sp library - spatial objects and methods in R
• GDAL libary - readOGR, writeOGR
• Convert Data Frame to SPDF
• CRS definitions and transformations
• Map Overlays
• Geospatial Maps with tmap